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Abstract 

In a network of dynamical systems, concurrent synchronization is a regime 
where multiple groups of fully synchronized elements coexist. In the brain, con- 
current synchronization may occur at several scales, with multiple "rhythms" in- 
teracting and functional assemblies combining neural oscillators of many different 
types. Mathematically, stable concurrent synchronization corresponds to conver- 
gence to a flow-invariant linear subspace of the global state space. We derive a 
general condition for such convergence to occur globally and exponentially. We 
also show that, under mild conditions, global convergence to a concurrently syn- 
chronized regime is preserved under basic system combinations such as negative 
feedback or hierarchies, so that stable concurrently synchronized aggregates of ar- 
bitrary size can be constructed. Robustnesss of stable concurrent synchronization 
to variations in individual dynamics is also quantified. Simple applications of these 
results to classical questions in systems neuroscience and robotics are discussed. 

1 Introduction 



Distributed synchronization phenomena are the subject of intense research. In the brain, 
such phenomena are known to occur at different scales, and are heavily studied at 
both the anatomical and computational levels. In particular, synchronization has been 
proposed as a general principle for temporal binding of multisensory data jl21 ^1 
EOl im 123 E2] 5 and as a mechanism for perceptual grouping 51j, neural computation 
miHnj and neural communication [211 [13 EHl l,40j . Similar mathematical models describe 
fish schooling or certain types of phase-transition in physics |45j . 
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In an ensemble of dynamical elements, concurrent synchronization is defined as a 
regime where the whole system is divided into multiple groups of fully synchronized 
elements^, but elements from different groups are not necessarily synchronized Pl l52[ IHB] 
and can be of entirely different dynamics ^T] . It can be easily shown that such a regime 
corresponds to a flow-invariant linear subspace of the global state space. Concurrent 
synchronization phenomena are likely pervasive in the brain, where multiple "rhythms" 
are known to coexist pT| IHH] , neurons can exhibit many qualitatively different types of 
oscillations fU] ITI)] , and functional models often combine multiple oscillatory dynamics. 

In this paper, we introduce a simple sufficient condition for a general dynamical 
system to converge to a flow-invariant subspace. Our analysis is built upon nonlinear 
contraction theory |221|iBj, and thus it inherits many of the theory's features : 

• global exponential convergence and stability are guaranteed, 

• convergence rates can be exphcitly computed as eigenvalues of well-defined sym- 
metric matrices, 

• robustness to variations in dynamics can be easily quantified, 

• under simple conditions, convergence to a concurrently synchronized state can be 
preserved through system combinations. 

As we shall see, under simple conditions on the coupling strengths, architectural 
symmetries and/or diffusion- like couplings create globally stable concurrent syn- 
chronization phenomena. This is illustrated in figure 1, which is very loosely inspired 
by oscillations in the thalamocortical system |2S1 112 ESI IS21 EH] ■ Qualitatively, global 
stability of the concurrent synchronization is in the same sense that an equilibrium point 
is globally stable - any initial conditions will lead back to it, in an exponential fashion. 
But of course it can yield extremely complex, coordinated behaviors. 

Section 2 recalls key concepts of nonlinear contraction theory and derives a theoreti- 
cal tool for studying global convergence to a flow-invariant subspace. Section 3 presents 
the paper's main mathematical results, relating stable concurrent synchronization to 
coupling structures and flow-invariant subspaces created by symmetries or diffusion-like 
couplings. Robustnesss of concurrent synchronization to variations in individual dy- 
namics is also quantified, showing in particular how approximate symmetries lead to 
quasi-synchronization. Section 4, motivated by evolution and development, studies con- 
ditions under which concurrent synchronization can be preserved through combinations 
of multiple concurrently synchronized regimes. Finally, section 5 discusses potential 
applications of these results to general questions in systems neuroscience and robotics. 

^In the literature, this phenomenon is often called poly-, or cluster or partial synchronization. How- 
ever, the last term can also designate a regime where the elements are not fully synchronized but behave 
coherently 



2 



Figure 1: An example of concurrent synchronization. Systems and connections of the same 
shape (and color) have identical dynamics (except black arrows, which represent arbitrary 
connections, and dashed arrows, which represent diffusive connections). This paper shows 
that under simple conditions on the coupling strengths, the group of (green) squares globally 
exponentially synchronizes (thus providing synchronized input to the outer elements), and so 
does the group of (yellow) diamonds, regardless of the specific dynamics, connections, or inputs 
of the other systems. 

2 Basic Tools 

2.1 Nonlinear contraction theory 

This section reviews basic results of nonlinear contraction theory |27 ^ I ^ Hn | l48j. which is 
the main stability analysis tool used in the paper. Essentially, a nonlinear time-varying 
dynamic system will be called contracting if initial conditions or temporary disturbances 
are forgotten exponentially fast, i.e., if trajectories of the perturbed system return to 
their nominal behavior with an exponential convergence rate. It turns out that relatively 
simple algebraic conditions can be given for this stability-like property to be verified, 
and that this property is preserved through basic system combinations. 

While we shall derive global properties of nonlinear systems, many of our results 
can be expressed in terms of eigenvalues of symmetric matrices ^3]. Given a square 
matrix A, the symmetric part of A is denoted by A^. The smallest and largest eigen- 
values of As are denoted by Amin(A) and Amax(A). Given these notations, the ma- 
trix A is positive definite (denoted A > 0) if Amin(A) > 0, and it is negative definite 
(denoted A < 0) if Amax(A) < 0. Finally, a square matrix A(x, t) is uniformly posi- 
tive definite if 3/? > 0, Vx, Vt : Amin(A(x, t)) > /5, and it is uniformly negative definite if 
3/3>0,Vx, Vt:A^i,(A(x,t))<-/5. 

The basic theorem of contraction analysis, derived in j2Z], can be stated as: 

Theorem 1 (Contraction) Consider, in R", the deterministic system 

x = f(x,t) (1) 
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where f is a smooth nonlinear function. Denote the Jacobian matrix off with respect to 
its first variable by If there exists a square matrix 0(x, t) such that 0(x, t)'''0(x, t) 
is uniformly positive definite and the matrix 

is uniformly negative definite, then all system trajectories converge exponentially to a 
single trajectory, with convergence rate | sup^t Amax(F)| > 0. The system is said to be 
contracting, F is called its generalized Jacobian, and ©(x, t)^©(x, t) its contraction 
metric. 

It can be shown conversely that the existence of a uniformly positive definite metric 
M(x, t) = 0(x, t)'''0(x, t) with respect to which the system is contracting is also a 
necessary condition for global exponential convergence of trajectories ^7\. Furthermore, 
all transformations corresponding to the same M lead to the same eigenvalues for the 
symmetric part F^ of F jl^l, and thus to the same contraction rate | sup^f Ainax(F)|. 

In the linear time-invariant system is globally contracting if and only if it is 

strictly stable, and F can be chosen as a normal Jordan form of the system with the 
coordinate transformation to that form ^ . Contraction analysis can also be derived 
for discrete-time systems and for classes of hybrid systems |2H] • 

Finally, it can be shown that contraction is preserved through basic system combina- 
tions (such as parallel combinations, hierarchies, and certain types of negative feedback, 
see 123 for details), a property which we shall extend to the synchronization context in 
this paper (section 4). 

Theorem 2 (Contraction and robustness) Consider a contracting systemic = i{'x,t) , 
with = 1 and contraction rate A. Let Pi{t) be a trajectory of the system, and let P2{t) 
be a trajectory of the disturbed system 

X = f(x,t) + d(x,t) 

Then the distance R(t) between Pi{t) and P2{t) verifies R{t) < sup^t ||d(x, t)||/A after 
exponential transients of rate A. 

For a proof and generalisation of this theorem, see section 3.7 in [T7\ . 
2.2 Convergence to a flow-invariant subspace 

We now derive a simple tool upon which the analyses of this paper will be based. The 
derivation is inspired by the idea of "partial" contraction, introduced in , which con- 
sists in applying contraction tools to virtual auxiliary systems so as to address questions 
more general than trajectory convergence. 
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Consider again, in M", the deterministic system 

x = f(x,t) (2) 

where f is a smooth nonhnear function. Assume that there exists a flow-invariant linear 
subspace M. (i.e. a hnear subspace 7V1 such that Vt : i{M.,t) C M), which imphes 
that any trajectory starting in M. remains in M.. Let p = dim(AI), and consider an 
orthonormal basis (ei, . . . , e„) where the first p vectors form a basis of A4 and the last 
n — p a. basis of A^"*". Define an (n — p) x n matrix V whose rows are ej^_^, . . . , e,j[. V 
may be regarded as a projection^ on M.-^^ and it verifies ^^EOI : 

V^V + U^U = In vv^ = I„_j, ^eM ^ Vx = 

where U is the matrix formed by the first p vectors. 

Theorem 3 Consider a linear flow-invariant subspace M. and the associated orthonor- 
mal projection matrix V. A particular solution ^p(t) of system ^ converges exponen- 
tially to M. if the system 

y = Vf(VTy + uTUxp(t),t) (3) 

is contracting with respect to y. 

If the above contraction condition is fullfilled for all^p, then starting from any initial 
conditions, all trajectories of system (0j will exponentially converge to Ai. If furthermore 
all the contraction rates for ^ are lower-bounded by some A > 0, uniformly in Xp and 
in a common metric, then the convergence to M. will be exponential with rate A (see 
figure{^. 

Proof : Let Zp = Vxp. By construction, Xp converges to the subspace Ml if and 
only if Zp converges to 0. Multiplying (j21) by V on the left, we get 

Zp = Vf (V^Zp + U^Uxp, t) (4) 

From (0}, y(t) = Zp{t) is a particular solution of system (jHl). In addition, since 
U'''Uxp G Ai and the linear subspace is flow- invariant, one has f(U^Uxp) G = 
Null(V), and hence y{t) = is another particular solution of system Q. If system Q 
is contracting with respect to y, then all its solutions converge exponentially to a single 
trajectory, which implies in particular that Zp{t) converges exponentially to 0. 

The remainder of the theorem is immediate. □ 

Corollary 1 A simple sufficient condition for global exponential convergence to M. is 
that 

Of 

V— < uniformly (5) 

^For simplicity we shall call V a "projection" , although the actual projection matrix is in fact V^V. 
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Figure 2: Convergence to a linear flow-invariant subspace 



or more generally, that there exists a constant invertible transform on Ai^ such that 



df 

eV— V^e^^ < uniformly 
ax 



(6) 



Proof : The Jacobian of Q with respect to y is 



V 



-(vTy + UTUx,(t),t) 



V 



so that the result is immediate by applying Theorem ^ □ 
Remarks 

• Non-orthonormal bases. In practice, the subspace Ai is often defined by the 
conjunction of (n—p) linear constraints. In a synchronization context, for instance, 
each of the constraints may be, e.g., of the form x, = Xj where Xj and Xj are 
subvectors of the state x. This provides directly a (generally not orthonormal) 



basis (e(,+i,...,e'„) 
and which verifies 



of Ai , and thus a matrix V whose rows are e' 



■p+i 



V = TV 



with T an invertible {n — p) x [n — p) matrix. We have x G 
and 

.9f_^ _ ,^,df. 



ax 



< 



ax 



< 



V'x = 



(7) 



Consider for instance three systems of dimension m and two systems of dimension 
p, and assume that Ai = {xi = X2,X5 = — IOX4} is the synchronization subspace 
of interest (with Xj denoting the state of each individual system). One has directly 



V 



I 



m 





Ir 



■■■ 

lOL 





In 



Note however that the equivalence in equation ((7j) does not yield the same upper 
bound for the eigenvalues of the two matrices. Thus, in order to compute explicitly 
the convergence rate to Ai, one has to revert to the orthonormal version, using 
e.g. a Gram-Schmidt procedure [13] on the rows of V. 
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• More general invariant subspaces. This theorem can be extended straightforwardly 
to time- varying affine invariant subspaces of the form m(t)+A^ (apply the theorem 
to x(t) = x(t) — m{t)). Preliminary results have also been obtained for nonlinear 
invariant manifolds 



2.3 Global synchronization in networks of coupled identical dy- 
namical elements 

In this section, we provide by using theorem El a unifying and systematic view on several 
prior results in the study of synchronization phenomena (see e.g. [ini EEl ESI ESI EE] and 
references therein). 

Consider first a network containing n identical dynamical elements with diffusive 
couplings (4S| 

±i = f (xj, t) + ^ Kij-(xj - Xi) i = l,...,n (8) 

Let L be the Laplacian matrix of the network (La = ^j-^jKjj, Lj^ = — Kjj for 
i ^ i), and'^ 



/xi \ 



f(x,t) 



Equation (jSj) can be rewritten in matrix form 



/ f(xi,t) \ 
\ f(x„,t) / 



X = f (x, t) — Lx 
The Jacobian matrix of this system is J = G — L, where 



(9) 



/ £(xi,t) 



G(x,t) 

















(hi 



Let now (ei, . . . , e^) be a basis of the state space of one element and consider the 
following vectors of the global state space 











ei = 




, • • • , Srf = 






[e, ) 







Let A4 = span{ei, . . . ,2^} be the "diagonal" subspace spanned by the e^. Note that 
X* e if and only if x^ = . . . = x* , i.e. all elements are in synchrony. In such a 

^The overscript ^ denotes a vector in the global state space, obtained by grouping together the 
states of the elements. 
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case, all coupling forces equal zero, and the individual dynamics are the same for every 
element. Hence 

/ f(x*i,t) \ 



f(S*,t) -Lx* 



G M 



\ f(x*i,t) J 



which means that A4 is flow-invariant. 

Consider, as in section 1?!^ the projection matrix V on A4^. Since V is built from or- 
thonormal vectors, Ainax(VG(x, t)V^) is upper-bounded by maxj Amax (^(xj, t)) . Thus, 
by virtue of theorem El a simple sufficient condition for global exponential synchroniza- 
tion is 

A^i„(VLV^) > supA^ax(^(a,t)^ (10) 

Furthermore, the synchronization rate, i.e. the rate of convergence to the synchro- 
nization subspace, is the contraction rate of the auxiliary system Q. 
Let us now make some brief remarks. 

(i) Undirected'^ diffusive networks. In this case, it is well known that L is symmetric 
positive semi-definite, and that is a subset of the eigenspace corresponding to 
the eigenvalue [Zj. Furthermore, if the network is connected, this eigenspace is 
exactly A^, and therefore VLV^ is positive definite (its smallest eigenvalue is called 
the network's algebraic connectivity 0). Assume now that L is parameterized by 
a positive scalar k (i.e. L = fcLo, for some Lq), and that ^ is upper-bounded. 
Then, for large enough k (i.e. for strong enough coupling strength), all elements 
will synchronize exponentially. 

(ii) Network of contracting elements. If the elements x, are already contracting when 
taken in isolation (i.e. § is uniformly negative definite), then in presence of weak 
or non-existent couplings (VLV^ = 0), the Jacobian matrix J of the global system 
will remain uniformly negative definite jlH]. Thus, the projected Jacobian matrix 
will be a fortiori uniformly negative definite, implying exponential convergence to 
the synchronized state. 

One can also obtain this conclusion by using a "pure" contraction analysis. Indeed, 
choose a particular initial state where xi(0) = . . . = x„(0). The trajectory starting 
with that initial state verifies Vt,xi(t) = . . . = x„(t) by fiow-invariance. Since the 
global system is contracting, any other initial conditions will lead exponentially 
to that particular trajectory, i.e., starting with any initial conditions, the system 
will exponentially converge to a synchronized state. 
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•'Undirected" is to be understood here in the graph-theoretical sense, i.e. : for all i,j, the connection 
from i to j is the same as the one from j to i. Therefore, an undirected network can be represented by 
an undirected graph, where each edge stands for two connections, one in each direction. 
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(iii) Nonlinear couplings. Similarly to (4^, the above result actually extends to non- 
linear couplings described by a Laplacian matrix L(x, t). Replacing the auxiliary 
system by 



(iv) Leader- followers network. Assume that there exists a leader in the network |48j . 
i.e., an element which has no incoming connections from the other elements, = 
f(x^,t). Convergence to ^A (guaranteed by satisfying (HHI)) then implies that all 
the network elements will synchronize to the leader trajectory x^(t). 

(v) Non-diffusive couplings. Note that the above results are actually not limited to 
diffusive couplings but apply to any system of the general form Q. This point 
will be further illustrated in sections 13.11 and 13.21 

3 Main discussion 

In nonlinear contraction theory, the analysis of dynamical systems is greatly simplified 
by studying stability and nominal motion separately. We propose a similar point of view 
for analyzing synchronization in networks of dynamical systems. In section ITU we study 
specific conditions on the coupling structure which guarantee exponential convergence 
to a linear subspace. In section we examine how symmetries and/or diffusion-like 
couplings can give rise to specific fiow-invariant subspaces corresponding to concurrent 
synchronized states. 

3.1 Some coupling structures and conditions for exponential 
synchronization 

3.1.1 Balanced diffusive networks 

A balanced network |33^ is a directed diffusive network which verifies the following 
equality for each node i (see figure El for an example) 



Because of this property, the symmetric part of the Laplacian matrix of the net- 
work is itself the Laplacian matrix of the underlying undirected graph to the network^. 

""In fact, it is easy to see that the symmetric part of the Laplacian matrix of a directed graph is the 
Laplacian matrix of some undirected graph if and only if the directed graph is balanced. 



y = Vf (V^y + U^Uip, t) - VL(Sp, t) (V^y + U^USp) 
the same steps show that global synchronization is achieved exponentially for 
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Thus, the positive definiteness of VLV^ for a balanced network is equivalent to the 
connectedness of some well-defined undirected graph. 



Figure 3: A balanced network with Lapla- 
cian matrix 
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Its underlying undirected graph, with 
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For general directed diffusive networks, finding a simple condition implying the pos- 
itive definiteness of VLV^ (such as the connectivity condition in the case of undirected 
networks) still remains an open problem. However, given a particular example, one can 
compute VLV^ and determine directly whether it is positive definite. 

3.1.2 Extension of diffusive connections 

In some applications 49 , one might encounter the following dynamics 

Xl = fi(xi,t) + A;A'^(Bx2 - Axi) 
±2 = f2(x2,t) + A;BT(Axi - BX2) 

Here xi and X2 can be of different dimensions, say di and d2- A and B are constant 
matrices of appropriate dimensions. The Jacobian matrix of the overall system is 



dfi 
(9xi 



9x2 



— /cL, where L 



A^A -A^B 
-B^A B^B 



Note that L is symmetric positive semi-definite. Indeed, one immediately verifies 



that 



Vxi,X2 : ( Xl X2 ) L ^ = (Axi - Bx2)^(Axi - BX2) > 

Consider now the linear subspace of M^^ x M'^^ defined by 

= I ( I G M'^i X M'^^ . Axi - Bx2 = 
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and use as before the orthonormal projection V on Ai^, so that VLV^ is positive 
definite. Assume furthermore that M is flow-invariant, i.e. 

V(xi,X2) G M^^ X R'^^ [Axi = Bx2] ^ [Afi(xi) = Bf2(x2)] 

and that the Jacobian matrices of the individual dynamics are upper-bounded. Then 
large enough k, i.e. for example 

fcAmin(VLV^) > max sup Amax^(ai, t) 

ensures exponential convergence to the subspace A4. 

The state corresponding to Ai can be viewed as an extension of synchronization 
states to systems of different dimensions. Indeed, in the case where xi and X2 have the 
same dimension and where A = B are non singular, we are in the presence of classical 
diffusive connections, which leads us back to the discussion of section IT^ 

As in the case of diffusive connections, one can consider networks of so-connected 
elements, for example : 

fi(xi,t) + AJ(BaX2 - Aijxi) + AJ(CaX3 - Ac-xi) 

f2(x2, t) + Bj(CijX3 - BcX2) + B^(AbXi - B^X2) 
f3(x2,t) + CT(AcXi - C^Xs) + C5(BcX2 - CbXs) 

leads to a positive semi-deflnite Laplacian matrix 

\ / AlAc -AJCa 

BjBc -BJCb + 

-C]jBc CJCb J \ —C\Ac C\Ca 

and potentially a flow-invariant subspace 

M = {Asxi = B^X2} n {BcX2 = Cijxa} n {C^xg = AcXi} 

The above coupling structures can be implemented in nonlinear versions of the pre- 
dictive hierarchies used in image processing (e.g. 1^ l?f| 1^ ITUl 1^ ). 

3.1.3 Excitatory-only networks 

One can also address the case of networks with excitatory-only connections. Consider 
for instance the following system and its Jacobian matrix^ 

xi = fixi,t) + kx2 T ff(a;i,t) 1 
i2 = /(x2,t) + A;xi 9l^X2,t)J^ \1 




AIAb 


-A^B^ 








-B\Ab 


B\Ba 
















I 




^ 



^For the sake of clarity, the elements are assumed to be 1-dimensional. However, the same reasoning 
applies for the multidimensional case as well: instead of span| ( i ) |i o^^s considers spanjei, . . . , e^} 
as in section IT^ 
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Clearly, span{(l,l)} is flow-invariant. Applying the methodology described above, 
we choose V = -^(1, —1), so that the projected Jacobian matrix is | (ff (^^i, t) + ^{x2, t)) — 

k. Thus, for k > sup^_^ |^(a,t), the two elements synchronize exponentially. 

In the case of diffusive connections, once the elements are synchronized, the coupling 
terms disappear, so that each individual element exhibits its natural, uncoupled behav- 
ior. This is not the case with excitatory-only connections. This is illustrated in flgureE] 
using FitzHugh-Nagumo oscillator models (see appendix^ for the contraction analysis 
of coupled FitzHugh-Nagumo oscillators). 




Figure 4: From left to right : a single oscillator, two oscillators coupled through diffusive 
connections, two oscillators coupled through excitatory-only connections. 



3.1.4 Rate models for neuronal populations 

In computational neuroscience, one often uses the following simplifled equations to model 
the dynamics of neuronal populations 

T±i = -Xj + $ kijXj{t)^ + Ui{t) 

Assume that the external inputs Uj(t) are all equal, and that the synaptic con- 
nections kij verify 3c, Vi, "^j^ikij = c (i.e., that they induce input-equivalence, see 
section IT^ . Then the synchronization subspace {xi = . . . = x„} is flow- invariant. Fur- 
thermore, since each element, taken in isolation, is contracting with contraction rate 
1/r, synchronization should occur when the coupling is not too strong (see remark (ii) 
in section IT^ . 

Speciflcally, consider flrst the case where $ is a linear function : $(x) = yux. The 
Jacobian matrix of the global system is then — 1„ -|- /xK, where K is the matrix of kij. 
Using the result of remark (ii) in section Y2.'S\ a sufficient condition for the system to be 
contracting (and thus synchronizing) is that the couplings are weak enough (or more 
precisely, such that /iAmax(K) < 1). 

The same condition is obtained if $ is now e.g. a multidimentional sigmoid of 
maximum slope /i (see remark (iii) in section 1231) • 

Besides the synchronization behavior of these models, their natural contraction prop- 
erty for weak enough couplings of any sign is interesting in its own right. Indeed, given 
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a set of (not necessarily equal) external inputs Uj(t), all trajectories of the global system 
will converge to a unique trajectory, independently of initial conditions. 

3.2 Symmetries, diffusion-like couplings, flow-invariant sub- 
spaces and concurrent synchronization 

Synchronized states can be created in at least two ways : by architectural and internal'' 
symmetries HH [HI El EHI or by diffusion- like couplings gHl 1201 E3 123 12] • Actually, we 
shall see that both, together or separately, can create flow-invariant subspaces corre- 
sponding to concurrently synchronized states. 

3.2.1 Symmetries and input-equivalence 

In section 1221 we argued that, in the case of coupled identical elements, the global syn- 
chronization subspace Ai represents a flow-invariant linear subspace of the global state 
space. However, several previous works have pointed out that larger (less restrictive) 
flow-invariant subspaces may exist if the network exhibits symmetries [^21 1211201! even 
when the systems are not identical 

The main idea behind these works can be summarized as follows. Assume that the 
network is divided into k aspiring synchronized groups Si, . . . , S^^- The flow-invariant 
subspace corresponding to this regime (in the sequel, we shall call such a subspace a 
concurrent synchronization subspace), namely 

{(xi; . . . ;x„) : VI < m < fc,Vi, j G 5^ : Xj = x^} 

is flow- invariant if, for each Sm, the following conditions are true : 

(i) if i,j G Sm, then they have a same individual (uncoupled) dynamics 

(ii) if i, j G Sm, and if they receive their input from elements i' and j' respectively, then 
i' and j' must be in a same group Sm' , and the coupling functions (the synapses) 
i' —>■ i and j' — > j must be identical. If i and j have more than one input, they must 
have the same number of inputs, and the above conditions must be true for each 
input. In this case, we say that i and j are input-symmetric, or more precisely, 
input- equivalent (since formally "symmetry" implies the action of a group). 

One can see here that symmetry, or more generally input-equivalence, plays a key 
role in concurrent synchronization. For a more detailed discussion, the reader is referred 

to [II1II2]. 

^Internal symmetries can easily be analyzed within our framework as leading to flow-invariant sub- 
spaces, and we shall use this property in section IKl^ for building central pattern generators. However, 
they will not be discussed in detail in this article. The interested reader can consult [H]. 

^Some groups may contain a single element, see section 14.2.31 
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Remark : One can thus turn on/ofF a specific symmetry by turning on/off a single 
connection. This has similarities to the fact that a single inhibitory connection can turn 
on/off an entire network of synchronized identical oscillators |I5] . 

3.2.2 Diffusion-like couplings 

The condition of input-equivalence can be relaxed when some connections within a 
group are null when the connected elements are in the same state. Such connections are 
pervasive in the literature : diffusive connections (in a neuronal context, they correspond 
to electrical synapses mediated by gap junctions imE], in an automatic control context, 
they correspond to poursuit or velocity matching strategies [331 120]; ■ ■ ■), connections 
in the Kuramoto model (IHl 1201 HHI (i-e. in the form Xj = /(xj, t) + kij sm{xj — Xj)), 
etc. 

Indeed, consider for instance diffusive connections and assume that 

• i' ^ i has the form Ki(xj/ — Xj) 

• j' ~^ j has the form K2(xj/ — Xj) with possibly Ki ^ K2 

Here, i and j are not input-equivalent in the sense we defined above, but the subspace 
{xj = Xj = Xj/ = Xj/} is still flow-invariant. Indeed, once the system is on this synchro- 
nization subspace, we have Xj = Xj/, Xj = Xj/, so that the diffusive couplings i' ^ i and 
j' — >■ j vanish. 

One can also view the network as a directed graph G, where the elements are repre- 
sented by nodes, and connections i j hj directed arcs i ^ j. Then, the above remark 
can be reformulated as 

1 : for all m, color the nodes of Sm with a color m, 

2 : for all m, erase the arcs representing diffusion- like connections and joining two nodes 

in Syn, 

3 : check whether the initial coloring is balanced (in the sense of [TT|) with respect to 

the so-obtained graph. 

It should be clear by now that our framework is particularly suited to analyze con- 
current synchronization. Indeed, a general methodology to show global exponential 
convergence to a concurrent synchronization regime consists in the following two steps 

• First, find an flow-invariant linear subspace by taking advantage of potential sym- 
metries in the network and/or diffusion-like connections. 
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• Second, compute the projected Jacobian matrix on the orthogonal subspace and 
show that it is uniformly negative definite (by explicitly computing its eigenvalues 
or by using results regarding the form of the network, e.g. remark (i) in section 
12.31 or section EIH). 

3.3 Illustrative examples 




Figure 5: Three example networks 



(i) The first network has three non-trivial fiow-invariant subspaces other than the 
global sync subspace, namely A4i = {xi = X2,X3 = X4}, A^2 = {^1 = ^3,^2 = 
X4}, and A4s = {xi = X4,X2 = X3}. Any of these subspaces is a strict superset 
of the global sync subspace, and therefore one should expect that the convergence 
to any of the concurrent sync state is "easier" than the convergence to the global 
sync state [^S 121 EE]- This can be quantified from (fTO|) . by noticing that 

MA^Ms^M^CMi X^U'yA'LVl) > Xn.in(VBLYl) (11) 

While in the case of identical systems and relatively uniform topologies, this "per- 
colation" effect may often be too fast to observe, (fTT|) applies to the general concur- 
rent synchronization case and quantifies the associated and possibly very distinct 
time-scales. 

(ii) The second network has only one non-trivial fiow-invariant subspace {xi = X2, X3 = 
X4}. 

(iii) If the dashed blue arrows represent diffusive connections then the third network 
will have one non-trivial flow-invariant subspace {x2 = X3 = X4,X5 = xg = xy}, 
even if these extra diffusive connections obviously break the symmetry. 

Let's study in more detail this third network, in which the connections between the 
round element and the square ones are modelled by trigonometric functions (we shall 
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see in section I4.2.HI that their exact form has no actual influence on the convergence 
rate). 

' ' + 02 sin(t;3) 

+ CiVq 

+ 61 (f 2 - fs) 
+ 63(f3 -1-4) 

diVe) 
divr) 



V3 
V4 

V7 
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The Jacobian matrix of the couplings is 
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As we remarked previously, the concurrent synchronization regime {f 2 = f 3 = f 4, f 5 = 
^6 = ^7} is possible. Bases of the linear subspaces M. and M.-^ corresponding to this 
regime are 
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Group together the vectors of the basis of M. 
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into a matrix V and compute 
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As a numerical example, let h\ = 3a, 62 = 4a, 63 = 5a, ci = a, C2 = 2a, di = 3a, 
d2 = 4:a and evaluate the eigenvalues of VLgV^. We obtain approximately 1.0077a 
for the smallest eigenvalue. Using again FitzHugh-Nagumo oscillators and based on 
their contraction analysis in appendix ^ concurrent synchronization should occur for 
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a > 10.25. A simulation is shown in figure IHl One can see clearly that, after a transient 
period, oscillators 2, 3, 4 are in perfect sync, as well as oscillators 5, 6, 7, but that the 
two groups are not in sync with each other. 




Figure 6: Simulation result for network 3. 



3.4 Robustness of synchronization 

So far, we have been considering exact synchronization of identical elements. However 
this assumption may seem unrealistic, since real systems are never absolutely identical. 
We use here the robustness result for contracting systems (see theorem |2) to guarantee 
approximate synchronization even when the elements are not identical. 
Consider, as in section 12. 3[ a network of n dynamical elements 

Xj = fi(xi,t) + ^Kij(xj -Xi) i = l,...,n (12) 
with now possibly fj 7^ fj for i j. This can be rewritten as 
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c(x„,t) j 




x„ 
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-c(x„,t) j 



where c is some function to be defined later. Keeping the notations introduced in section 
12. 3[ one has 

X = c(x, t) — Lx + d(x, t) 

where d(x, t) stands for the last term of equation (fTSj) . 
Consider now the projected auxiliary system on 

y = Vc(VV + UU^S, t) - VLV^y + Vd(V^y + UU"^S, t) (14) 

Assume that the connections represented by L are strong enough (in the sense of 
equation (|Tn|l ). so that the undisturbed version of (|THl is contracting with rate A > 0. 
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Let D = supxf ||Vd(x, where D can be viewed as a measure of the dissimilarity of 
the elements. Since y = is a particular solution of the undisturbed system, theorem 
121 implies that the distance R{t) between any trajectory of (|14p and verifies, after 
a transient period, R{t) < D/X. In the x-space, it means that any trajectory will 
eventually be contained in a boundary layer of thickness D/X around the synchronization 
subspace M.. 

The choice of c can now be specified so as to minimize D/X. Neglecting for simplicity 
the variation of A, a possible choice for c(x, t) is then the center of the ball of smallest 
radius containing fi(x, t), . . . , f„(x, t), with D being the radius of that ball. 

Consider for instance, the following system (similar to the model used for coincidence 
detection in |1H] and section EH]) 

= f{xi) + k{xo - Xi) where I^in < h< ^max, Vz 

In this case, choosing c(x) = f{x) + -^max+Jmin ^ ^g^^ achieve the bound D/X, where A 
is the contraction rate of / and D = 

Remark : Assume that two spiking neurons are approximately synchronized, as 
just discussed. Then, since spiking induces large abrupt variations, the neurons must 
spike approximately at the same time. More specifically, if the bound on their trajectory 
discrepancy guaranteed by the above robustness result is significantly smaller than spike 
size, then this bound will automatically imply that the two neurons spike approximately 
at the same time. 

4 Combinations of concurrently synchronized groups 

This section shows that, under mild conditions, global convergence to a concurrently 
synchronized regime is preserved under basic system combinations, and thus that stable 
concurrently synchronized aggregates of arbitrary size can be systematically constructed. 
The results, derived for two groups, extend by recursion to arbitrary numbers of groups. 

4.1 The input-equivalence preservation assumption 

Consider two independent groups of dynamical elements, say 5*1 and 5*2. For each group 
i {i = 1,2), assume that a flow-invariant subspace Aii corresponding to a concurrently 
synchronized regime exists. Assume furthermore that contraction to this subspace can 
be shown, i.e. VjJjV^ < for some projection matrix Vj on A^^. 

Connect now the elements of 5*1 to the elements of 5*2, while preserving input- 
equivalence for each aspiring synchronized subgroup of Si and ^2. Thus, the combined 
concurrent synchronization subspace TVli x A^2 remains a flow-invariant subspace of the 

new global space. A projection matrix on (A^i x Ai2)^ can be V = 

where each Vj has been rescaled into Wj to preserve orthonormality. 



Wi 
W2 
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Specific mechanisms facilitating input-equivalence preservation will be discussed in 
section I4.2.HI 



4.2 Typology of combinations 

Let us now study several combination operations of concurrently synchronized groups 
and discuss how they can preserve convergence to a combined concurrent sync state. 

In 14.2.11 and I4.2.2[ the input-equivalence preservation condition of section 14.11 is im- 
plicitly assumed, and the results reflect similar combination properties of contracting 
systems j^ZlllSllini- More generally, as long as input-equivalence is preserved, any com- 
bination property for contracting systems can be easily "translated" into a combination 
property for synchronizing systems. 



4.2.1 Negative feedback combination 

The Jacobian matrices of the couplings are of the form J12 = —kjj^, with k a positive 
constant. Thus, the Jacobian matrix of the global system can be written as 




As in equation ^ of section 1221 consider a transform over {A4i x ^^2)"*" 

The corresponding generalized projected Jacobian matrix on {Aii x 7^2)^ is 
e VJV^ = { ^ ^ ^ ^ ^ < uniformly 

^ ^ \^ V^W2J2lWiT W2J2W2^ J 

so that global exponential convergence to the combined concurrent synchronization state 
can be then concluded. 



4.2.2 Hierarchical combination 

Assume that the elements in Si provide feedforward to elements in 5*2 but do not 
receive any feedback from them. Thus, the Jacobian matrix of the global system is 
Ji 



J21 J2 



. Assume now that W2 J21W1 is uniformly bounded and consider the 

I 

el 



coordinate transform = ( ^ ^ ) • Compute the generalized projected Jacobian 



\ eW2J2iWi W2J2W2 / 
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Since W2J21W1 is bounded, and Wi(Ji)Wi and W2(J2)W2 are both negative 
definite, 0e(VJV^)07^ will be negative definite for small enough e. 

Note that classical graph algorithms |22] allow large system aggregates to be system- 
atically decomposed into hierarchies (directed acyclic graphs, feedforward networks) of 
simpler subsystems (strongly connected components) [IH]. Input-equivalence then needs 
only be verified top-down. 

4.2.3 Case where Si has a single element 

Denote this element by ei (figure Q shows such a configuration where ei is the round red 
central element, and where 5*2 is the set of the remaining elements). Connections from 
(resp. to) ei will be called 1^2 (resp. 2— >1) connections. Then some simplifications 
can be made : 

• Input-equivalence is preserved whenever, for each aspiring synchronized subgroup 
of 5*2, the l—>-2 connections are identical for each element of this subgroup (in 
figure n the connections from ei to the yellow diamond elements are the same). 
In particular, one can add/suppress/ modify any 2—>-l connection without altering 
input-equivalence . 

• Since dim(A1]'-) = (a single element is always synchronized with itself), one 
has [Ail X Ai2)^ = ■^2- Thus, concurrent synchronization (and the rate of 
convergence) of the combined system only depends on the parameters and the 
states of the elements within 5*2. In particular, it neither depends on the actual 
state of Ci, nor on the connections from/to Ci (in figure^ the black arrows towards 
the red central element are arbitrary). 

In practice, the condition of identical 1— s>2 connections is quite pervasive. In a neu- 
ronal context, one neuron with high fan-out may have lO'^ identical outgoing connections. 

It is therefore quite easy to preserve an existing concurrent synchonization behavior 
while adding groups consisting of a single element. Thus, stable concurrent synchro- 
nization can be easily built one element at a time. 

4.2.4 The diffusive case 

We stick with the case where 5*1 consists of a single element, but make now the addi- 
tionnal requirement that the l—>-2 connections and the internal connections within 5*2 
are diffusive (so far in this section |3J we have implicitly assumed that the connections 
from Si to Sj only involve the states of elements in Si). The Jacobian matrix of the 
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combined system is now of the form 



/ dg{xi,t) 
dxi 



\ 



dfi(x2,t) 

8X2 



dfq{Xn,t) 
dXn 



( * 



-A:i 
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kq 





Lint 



where the first matrix describes the internal dynamics of each element, the second, the 
diffusive connections between ei and 5*2 (where 5*2 has q aspiring synchronized sub- 
groups), and the third, the internal diffusive connections within 32- 
Hence, the projected Jacobian matrix on {Aii x 



/ dfl(x2,t} 

' 8x2 



\ 



\ 



8fq(Xn,t) 

8x„ 



V' - v. 



/ kl 



\ 



is 



kn 



An interpretation of this remark is that there are basically three ways to achieve 
concurrent synchronization within 5*2, regardless of the behavior of element Ci and of its 
connections : 



(i) one can increase the strengths ki,...,kq of the 1^2 connections (which corre- 
sponds to adding inhibitory damping to S2), so that each element of 5*2 becomes 
contracting. In this case, all these elements will synchronize because of their con- 
tracting property even without any direct coupling among them (Ljnt = 0) (this 
possibility of synchronization without direct coupling is exploited in the coinci- 
dence detection algorithm of jHO], and again in section IKTTl of this paper), 

(ii) or one can increase the strength Lint of the internal connections among the elements 
of ^2, 

(iii) or one can combine the two. 



4.2.5 Parallel combination 

The elementary fact that, if VJjV^ < for a set of subsystem Jacobian matrices Jj, 
then V (^j ai(t) Jj) < for positive ai(t), can be used in many ways. Note that 
it does not represent a combination of different groups as in the above paragraphs, but 
rather a superposition of different dynamics within one group. 

One such interpretation, as in for contracting systems, is to assume that for a 
given system x = f(x, t), several types of additive couplings Lj(x, t) lead stably to the 
same invariant set, but to different synchronized behaviors. Then any convex combina- 
tion {ai(t) > 0, J2i (^i(t) = 1) of the couplings will lead stably to the same invariant set. 
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Indeed, 

f (x, t)-J2 WLi(x, t) = J2 W [f ^) - Li(x, t)] < 

The Lj(x, t) can be viewed as synchronization primitives to shape the behavior of 
the combination. 

5 Examples 

In conclusion, let us briefly discuss some general directions of application of the above 
results to a few classical problems in systems neuroscience and robotics. 

5.1 Coincidence detection for multiple groups 

Coincidence detection is a classic mechanism proposed for segmentation and classifica- 
tion. In an image for instance, elements moving at a common velocity are typically 
interpreted as being part of a single object, and this even when the image is only com- 
posed of random dots [23 12^1 ■ 

As mentioned in section 13^21 the possibility of decentralized synchronization via cen- 
tral diffusive couplings can be used in building a coincidence detector. In inspired 
in part by P], the authors consider a leader-followers network of FitzHugh-Nagumo^ 
oscillators, where each follower oscillator i (an element of 5*2, see section l4. 2. 4j) receives 
an external input J, as well as a diffusive coupling from the leader oscillator (the element 
ei of Si). Oscillators i and j receiving the same input (/, = Ij) synchronize, so that 
choosing the system output as ^i<j<„['i'i]''' captures the moment when a large number 
of oscillators receive the same input. 

However, the previous development also implies that this very network can detect 
the moments when several groups of identical inputs exist. Furthermore, it is possible 
to identify the number of such groups and their relative size. Indeed, assume that the 
inputs are divided into k groups, such that for each group Sm, one has Vz, j G Sm, h = Ij- 
Since the oscillators in Sm only receive as input (a) the output of the leader, which is the 
same for everybody and (b) the external input Jj, which is the same for every oscillator 
in group Sm, they are input-symmetric and should synchronize with each other (cf. 
section ISI21 and section E. 2. 4^ . 

Some simulation results are shown in figure [7| Note that contrast between groups 
could be further enhanced by using nonlinear "synapses", e.g. introducing input- 
dependent delays, which would preserve the symmetries. Similarly, any feedback mech- 
anism to the leader oscillator would also preserve the input-symmetries. 

Finally, adding all-to-all identical connections between the follower oscillators would 
preserve the input-symmetries and further increase the convergence rate, but at the 
price of vastly increased complexity. 

^See appendix 1X1 
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Figure 7: Simulation results for the coincidence detector. The network is composed of one 
leader with a constant input Iq = 25, and 80 followers whose inputs vary with time as shown 
in the upper figure. The lower figure plots the system output ^X^i<j<8o[^j]^) against time. 
One can clearly observe the existence of two successive, well separated spikes per period in a 
time interval around t = 250. Furthermore, one spike is about twice as large as the other one. 
This agrees with the inputs, since around t = 250, they are divided into two groups : 1/3 of 
them with value 37, 2/3 with value 30. 



5.2 Fast symmetry detection 

Symmetry, in particular bilateral symmetry, has also been shown to play a key role 
in human perception Consider a group of oscillators having the same individual 
dynamics and connected together in a symmetric manner. If we present to the network 
an input having the same symmetry, some of the oscillators will synchronize as predicted 
by the theoretical results of section 13.21 

One application of this idea is to build a fast bilateral symmetry detector (figures 
IHl El CHI), extending the oscillator-based coincidence detectors of the previous section. 
Although based on a radically different mechanism, this symmetry detector is also some- 
what reminiscent of the device in [Sj. 
Some variations are possible : 
(i) Other types of invariance. It is easy to modify the network in order to deal with 
multi-order (as opposed to bilateral) symmetry, or other types of invariance (trans- 
lation, rotation, . . . ). In each case, the network should have the same invariance 
pattern as what it is supposed to detect. 
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Figure 8: A fast bilateral symmetry detector. 
Every oscillator has the same dynamics as its mirror image and 
the two are coupled through diffusive connection. The inputs are 
provided to the detector through the bottom (red) layer. Then 
"information" travels bottom-up : each layer is connected to the 
layer right above it. Top-down feedback is also possible. 
Assume now that a mirror symmetric image is submitted to the 
network. The network, which is mirror symmetric by construc- 
tion, now receives a mirror symmetric input. Thus, the concur- 
rent synchronization subspace where each oscillator is exactly in 
the same state as its mirror image oscillator is flow-invariant. Fur- 
thermore, the diffusive connections, if they are strong enough (see 
2.2), guarantee contraction on the orthogonal space. By using the 
theoretical results above, one can deduce the exponential conver- 
gence to the concurrent synchronization regime. In particular, the 
difference between the top two oscillators should converge expo- 
nentially to zero. 

(ii) Since the exponential convergence rate is known, the network may be used to track 
time-varying inputs, as in the coincidence detection algorithm of |5Uj . 

(iii) Multidimensional inputs. Coincidence detectors and symmetry detectors may also 
handle multidimensional inputs. Two approaches are possible. One can either 
"hash" each multidimensional input into a one- dimensional input, and give the 
set of so-obtained one-dimensional inputs to the network. Or one can process each 
dimension independently in separate networks and then combine the results in a 
second step. 

5.3 Central pattern generators 

In an animal/robotics locomotion context, central pattern generators are often modelled 
as coupled nonlinear oscillators delivering phase-locked signals. We consider here a 
system of three coupled 2-dimensional Andronov-Hopf oscillators very similar to 
the ones used in the simulation of salamander locomotion ^H] : 

Xi = f (xi) + A;(R27rX2 — Xi) 
X2 = f (X2) + A;(R^X3 - X2) 
X3 = f (xg) + A;(R22rXi - X3) 
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Figure 9: Simulation on an artificial image. 

We create a 56 x 60 pixels symmetric image from a real picture 
of one of the authors. We give it as input to a network similar 
to the one in figure |H1 The first (bottom) layer of the network is 
composed of 7 x 6 = 42 FitzHugh-Nagumo oscillators (21 pairs) 
each receiving the sum of the intensities of 8 x 8 = 64 pixels, thus 
covering at every instant an active window of 56 x 48 pixels. The 
second layer consists of 4 oscillators, each receiving inputs from 9 
or 12 oscillators of the first layer. The third layer is composed of 
2 oscillators. 

At t = 0, the active window is placed on the left of the image (red 
box) and, as t increases, it slides towards the right. At t = T/2, 
where T is the total time of the simulation, the position of the 
window is exactly at the center of the image (green box) (see the 
simulation results in figure EJ. 



where f is the dynamics of an Andronov-Hopf oscillator and the matrix R22; describes 



a ^ planar rotation : 



X — y — X — xy 
X + y — — yx"^ 
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We can rewrite the dynamics as x = f(x) — A;Lx, where 
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First, observe that the linear subspace M. = | (^Rl^- (x), Ra^r (x), xj : x G IK' j- is 

fiow-invariant^°, and that A4 is also a subset of Null(Ls). Next, remark that the char- 
acteristic polynomial of is (X — 3/2)"^ so that the eigenvalues of are 0, with 
multiplicity 2, and 3/2, with multiplicity 4. Now since A4 is 2-dimensional, it is exactly 
the nullspace of L.,, which implies in turn that A^"*" is the eigenspace corresponding to 
the eigenvalue 3/2. 

Moreover, the eigenvalues of the symmetric part of ^ix,y) are 1 — (x^ + y'^) and 
1 — 3(x^ + y"^), which are upper-bounded by 1. Thus, for k > 2/3 (see equation (fTUj) in 
section IT^ . the three systems will globally exponentially converge to a ±^ -phase-locked 
state (i. e. a state in which the difference of the phases of two consecutive elements is 
constant and equals ±^). A computer simulation is presented in figure HH 

^°As it is suggested in footnote 7, the flow-invariance of A4 can be understood here as being "created" 
by the internal symmetries of the oscillators' dynamics. 
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Figure 10: Result for the simulation of figure IHl The figure shows \vi — V2\ where vi and 
V2 are the voltages of the two FN oscillators in the top layer. One can clearly observe that, 
around t = T/2 (T = 520 in this simulation), there is a short time interval during which the 
two oscillators are fully synchronized. 




Figure 11: Simulation for three coupled Andronov-Hopf oscillators. In the first two figures, we 
plot Ui against Xj for 1 < i < 3. Figure a shows the behavior of the oscillators for < t < 3s, 
figure b for 5.6s < t < 6s. In figures c and d, we plot xi, X2, X3 against time. Figure c for 
0<t< 0.5s, figure d for < t < 6s. 

Relaxing the symmetry or the diffusivity condition : In the previous example, 
the flow-invariance of the phase-locked state is due to (a) the internal symmetry of 
the individual dynamics f, (b) the global symmetry of the connections and (c) the 
"diffusivity" of the connections (of the form fc(Rx2 — xi)). Observe now, as in section 
13. 2[ that this flow-invariance can be preserved when one out of the two conditions (b) 
and (c) is relaxed. Consider for example the two following systems : 

• Symmetric but not "diffusive" : 

Xi = ffxi) + A;R27rX2 

3 

X2 = f(x2) + A;R^X3 

X3 = ffxs) + /cRa^rXi 

3 

(the connections are "excitatory-only" in the sense of section l!^.1.3j) . 
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• "Diffusive" but not symmetric 




f(xi) + A;i(RiX2 - Xl) 

f(x2) + /i;2(R2X3 - X2) 

f(x3) + A;3(R3Xi - X3) 



where the Rj represent any planar rotations such that R1R2R3 = I2 (i-e., any 
arbitrary phase-locking). 

By k;eeping in mind that for any planar rotation R and state x, one has f(Rx) = 
R(f (x)), it is immediate to show the flow-invariance of < ( Rl^. (x), Rg^r (x), x) : x G [ 



in the first case, and of {(RiR2(x), Ri(x), x) : x e M?} in the second case. Note however 
that the computations of the projected Jacobian matrices are different, and that in the 
first case the limit cycle's radius varies with k (cf. section 13. 1.3|) . 
Finally, note that 

(i) All the results of this section can be immediately extended to systems with more 



(ii) As compared to results based only on phase oscillators, this analysis guarantees 
global exponential convergence, rather than assuming that synchronization is al- 
ready essentially achieved. In addition, it exhibits none of the topological difficul- 
ties that may arise when coupling large numbers of phase oscillators. 

(iii) If f is less symmetric, only connections that exhibit the same symmetry as f can 
lead to a non-trivial flow-invariance subspace. 

(iv) It is also possible to extend this study to systems composed of oscillators with 
larger dimensions (living in for example), although a locomotion interpretation 
may be less relevant. 

5.4 Filtered connections and automatic gait selection 

Replacing ordinary connections in the CPG described in section 15. HI by filters enables 
frequency- based symmetry selection. This idea may have powerful applications, one of 
which could be automatic gait selection in locomotion. 

Consider for example the mechanism described in figure IT^ At low frequencies, the 
1 4-^ 3 and 2 <-> 4 connections are filtered out, so that the actual connections are 1^2 
and 3^4 (anti-synchronization) and 3^2 and 4 — 1 (quarter-period delay). The 
only non-trivial flow-invariant subspace is then {xi = R|(x3) = — X2 = R37r(x4)}. On 
the contrary, the 3^2 and 4^1 connections are filtered out at high frequencies, so 
that the flow-invariant subspace becomes {xi = X3 = — X2 = — X4}. 

Similarly to section 15.31 strong enough coupling gains ensure convergence to either 
of these two subspaces, according to the frequency at which the oscillators are running. 




3 



oscillators. 
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Figure 12: A CPG with filtered connections. 



The connections from the command box set the same 
frequency for the four oscillators. The 1 <-> 2 and 
3^4 arrows represent permanent anti- synchronization 
connections (i.e. connection j — > i is of the form 
k{—Xj — Xj)). The 1 <-> 3 and 2 <-> 4 arrows represent 
synchronization connections and they are high-pass fil- 
tered. Finally, the 3 — > 2 and 4^1 arrows stand for 
quarter-period delay connections (i.e. connection j ^ i 
is of the form k{Ii_^Xj — Xj), see section ESI and they 
are low-pass filtered. 



Note that standard techniques allow sharp causal filters with frequency-independent 
delays to be constructed easily 

An analogy with horse gaits could be made in this simplified setup, by associating 
the low- frequency regime with the walk (left fore, right hind, right fore, left hind), and 
the high-frequency regime with the trot (left fore and right hind simultaneously, then 
right fore and left hind simultaneously). Transitions between the two regimes would 
occur automatically according to the speed of the horse (the frequency of its gait). 

5.5 Temporal binding 

The previous development has suggested a mechanism for stable accumulation and in- 
teraction of concurrently synchronized groups, showing that the simple conditions for 
contraction to a linear subspace, combined with the high fan-out of typical neurons, 
increased the plausibility of large concurrently synchronized structures being created in 
the central nervous system in the course of evolution and development. The recently 
established pervasiveness of electrical synapses P would also be consistent with such 
architectures. 

More speculatively, different "rhythms" [a, j3, 7, 6) are known to coexist in the brain, 
which, in the light of the previous analysis, may be interpreted and modelled as concur- 
rently synchronized regimes. Since contracting systems driven by periodic inputs will 
have states of the same period j2ZI, different but synchronized computations could be 
robustly carried out by specialized areas in the brain using synchronized elements as 
their inputs. Such a temporal "binding" jl2llISll2SlllZII231E2linillllEl mechanism 
would also complement the general argument in that multisensory integration may 
occur through the interaction of contracting computational systems connected through 
an extensive network of feedback loops. In this context, and along the lines of section 
a translation to concurrent synchronization of recent results on centralized contracting 
combinations may be particularly relevant. Making these observations precise is the 
subject of future research. 
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A FitzHugh-Nagumo oscillators 

Some of our simulations involve coupled FitzHugh-Nagumo neural oscillators [HI |^ 

Vi = Vi{a - Vi){vi - 1) - Wi + li + k{vo - Vi) x < z < n 
Wi = (3vi --fWi - - 

In this paper, we use the following parameters values : a = 6, (3 = 3, 7 = 0.09. 
The contraction analysis of FitzHugh-Nagumo oscillators can be adapted from 
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